学习者&＃xff1a;骆栢维

题目来源&＃xff1a;生信基石之R语言

中级10 个题目&＃xff1a;http://www.bio-info-trainee.com/3750.html

备注:本文为笔者学习健明老师GitHub答案代码的学习收获&＃xff01;

作业1&＃xff1a;根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)

library(org.Hs.eg.db)symbol1&＃61;toTable(org.Hs.egSYMBOL)ensembl1&＃61;toTable(org.Hs.egENSEMBL)head(ensembl1)#查看前6行head(symbol1)#查看前6行exprset1 &＃39;gene_id&＃39;)#orexprset2 &＃39;gene_id&＃39;)c &＃39;ENSG00000000003.13&＃39;, &＃39;ENSG00000000005.5&＃39;, &＃39;ENSG00000000419.11&＃39;, &＃39;ENSG00000000460.15&＃39;, &＃39;ENSG00000000938.11&＃39;)c substr(c,c#看一下吧answer #这里按照逻辑值提取元素##来探讨一下吧exprset1$ensembl_id %in%c#%in%表示前者是否在后者里面&＃xff0c;返回逻辑值table(exprset1$ensembl_id %in%c)#ORc answer1 "ensembl_id",by.y&＃61;View(answer)

学习心得&＃xff1a;数据的提取

在上一篇分享文章的Q6&＃xff0c;笔者提及过&＃xff1a;提取数据的方法可以根据“名称”、“逻辑值” 和 “位置”。在上面的题目中&＃xff0c;笔者是通过逻辑值提取数据&＃xff0c;在下面的“读入数据” 中笔者也是采用逻辑值的方式提取数据。

作业2&＃xff1a;根据R包hgu133a.db找到下面探针对应的基因名(symbol)

install.packages("BiocManager")BiocManager::install(&＃39;hgu133a.db&＃39;)library(hgu133a.db)symbol2&＃61;toTable(hgu133aSYMBOL)#获取探针表格symbol2[1:4,]#看一下前4行c &＃39;1053_at&＃39;, &＃39;1320_at&＃39;, &＃39;1405_i_at&＃39;, &＃39;1431_at&＃39;, &＃39;1438_at&＃39; ,&＃39;1487_at&＃39; ,&＃39;1494_f_at&＃39;, &＃39;1598_g_at&＃39;, &＃39;160020_at&＃39;, &＃39;1729_at&＃39;,&＃39;177_at&＃39; )answer2 in% c,]table(symbol2$probe_id %in% c)#看一下结果#另一种方法同Q1

作业3&＃xff1a;找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量&＃xff0c;并且绘制在 progres.-stable分组的boxplot图

BiocManager::install(&＃39;CLL&＃39;)BiocManager::install(&＃39;hgu95av2.db&＃39;)library(CLL)data(package&＃61;"CLL")#查看R包里面具体数据集data("sCLLex")#加载R包内置数据,是列表biaoxing #提取表型数据exprSet #提取表达矩阵 library(hgu95av2.db)symbol3&＃61;toTable(hgu95av2SYMBOL)# 在symbol3中搜索TP53得到三个探针"1939_at"&＃xff1b;&＃39;31618_at&＃39;&＃xff1b;1974_s_at&＃39;probe in%probeexprSet boxplot(exprSet["1939_at",] ~ biaoxing$Disease,data &＃61; exprSet,col&＃61;"red") #其他的同上&＃xff0c;换基因名就可以了

作业3(补充)&＃xff1a;ggpubr

#ggpubr(以"1939_at"为例)library(ggpubr)#首先在绘图之前&＃xff0c;我们要处理数据&＃xff0c;把分组信息和表达信息合并起来#都是笔者发现他们的sampleID和表达数据的行名并不一致&＃xff0c;怎么办呢&＃xff1f;visual1 as.data.frame(exprSet)[biaoxing$SampleID2 function(x){ x paste(c(x,"CEL"),collapse &＃61; ".")#paste函数主要是为了将两个字符串合并起来}))##或者#############################################biaoxing$SampleID1 function(x){ x paste(c(x,"CEL"),collapse &＃61; ".")})##或者#######################################################上面是改分组信息名&＃xff0c;也可以改表达数据的列名exprSet as.data.frame(exprSet)colnames(exprSet) ".CEL",colnames(exprSet)############################################################identical(colnames(exprSet),biaoxing$SampleID)#判断表达数据列名和分组信息是否对应&＃xff0c;结果为Truebiaoxing$Disease as.character(biaoxing$Disease)b as.data.frame(t(biaoxing$Disease as.character(biaoxing$Disease)visual1 #按行合并visual1 % as.matrix()%>% t()%>% as.data.frame()visual1$&＃96;1939_at&＃96; as.numeric(ggboxplot(visual1, "Disease", "&＃96;1939_at&＃96;", color &＃61; "Disease", palette &＃61;c("#00AFBB", "#E7B800"), add &＃61; "jitter", shape &＃61; "Disease")

实战一下吧——lncRNA数据处理(模拟)

##########################实战一下吧############################a 1:b &＃39;tumor&＃39;, c "male",e ":")d "_")shizhan #创造一个数据框#1 191:normal_female#2 131:tumor_male#3 13:normal_male#4 266:normal_female#5 257:normal_female#6 246:normal_female###如果要把他们分三列怎么做呢library(tidyr)shizhan$d ":",shizhan "a",#完美############################################另一种方法shizhan1 "a",###发现宝藏&＃xff01;&＃xff01;&＃xff01;&＃xff01;&＃xff01;# a status gender#1 191 normal female#2 131 tumor male#3 13 normal male#4 266 normal female#5 257 normal female#6 246 normal female

学习心得

1.字符串的处理

*截取&＃xff1a;在上一篇分享文章中&＃xff0c;笔者介绍了unlist与strsplist截取字符串&＃xff0c;还有笔者也拓展了substr和substring的用法

*合并&＃xff1a;使用paste函数&＃xff0c;在这里要注意有时候collapse参数得不到理想的结果&＃xff0c;可以把参数变为sep&＃xff0c;二者都可以作为连接参数(具体为什么collapse有时候得不到理想的结果&＃xff0c;笔者也不确定)

*替换&＃xff1a;使用gsub函数&＃xff0c;可以做批量处理。

*拆分&＃xff1a;使用separate函数&＃xff0c;这里要注意的是&＃xff0c;它第一个参数为数据框&＃xff0c;第二个参数为列名&＃xff0c;不要觉得“数据框$列名” 就可以了。笔者就犯了这么愚蠢的错误。

2.apply ; lapply 和sapply的区别

*输入的数据&＃xff1a;apply为数据框和矩阵&＃xff1b;lapply为列表&＃xff0c;向量和数据框&＃xff1b;sapply同lapply。

*输出结果&＃xff1a;apply为向量&＃xff0c;列表和数组&＃xff1b;lapply为列表&＃xff1b;sapply为向量和矩阵。

*回到上面的合并字符串的前两种方法&＃xff0c;笔者在lapply前面加了unlist&＃xff0c;而sapply却没有。可以理解吗&＃xff1f;

温馨提示&＃xff1a;由于笔者的的网络不好&＃xff0c;没有办法下载所有数据(郁闷&＃xff01;&＃xff01;&＃xff01;)&＃xff0c;所以笔者按照作业相应的考点进行操作。接下来的内容可能与原作业存在差异&＃xff0c;还请大家海涵。

读入数据(对应作业8)

exprSet "GSE24673_series_matrix.txt",comment.char &＃61; library(tibble)rownames(exprSet) 1]exprSet -1]platformMap "platformMap.txt")platformMap[platformMap$gpl&＃61;&＃61;"GPL6244",3]#找到对应的注释包"hugene10sttranscriptcluster"if (!requireNamespace("BiocManager", quietly &＃61; TRUE)) install.packages("BiocManager")BiocManager::install("hugene10sttranscriptcluster.db")library(hugene10sttranscriptcluster.db)prode head(prode)

基因注释(对应作业4和9)

exprSet$probe_id #将行名转变为第一列&＃xff0c;笔者真是手残fexprset "probe_id")fexprset "probe_id",fexprset -1]fexprset 2:head(fexprset)#发现在第一列的是"RN7SK"基因maxgene "RN7SK",maxgenefexprset #去重rownames(fexprset) 1]final -1]save(final,file &＃61; "final_exprset.Rdata")#保存数据备注&＃xff1a;很感激果子老师为我们这些生信小白这里的平台文件的注释R包&＃xff01;详情请看&＃xff1a;https://mp.weixin.qq.com/s/nWbMO4mULgN__nPjooRDlg

差异分析(对应作业5、7和10一部分)

library(limma)group "tumor",group group,levels &＃61; c(design group)colnames(design) group)design# tumor normal#1 1 0#2 1 0#3 1 0#4 1 0#5 1 0#6 1 0#7 1 0#8 1 0#9 1 0#10 1 1#11 1 1fit #线性模型拟合fit2 #贝叶斯检验allDiff&＃61;topTable(fit2,adjust&＃61;&＃39;fdr&＃39;,coef&＃61;2,number&＃61;Inf) #输出结果save(allDiff,file &＃61; "allDiff.Rda")#保存dif_subset 1 & adj.P.Val % filter(adj.P.Val 1) %>% column_to_rownames()save(diffLab,group,file &＃61; "diffLab.Rda")##########################################################

热图绘制(对应作业6、7和10)

library(pheatmap)heatdata #筛选出差异表达明显的基因表达谱pheatmap(heatdata)#可以简单画一下&＃xff0c;都是没有分组信息group "tumor",annotation_col group)rownames(annotation_col) #设置分组信息pheatmap(heatdata, #热图的数据 annotation_col &＃61;annotation_col, #标注样本分类 annotation_legend&＃61;TRUE, # 显示注释 show_rownames &＃61; F,# 显示行名 scale &＃61; "row", #以行来标准化&＃xff0c;使结果更加明显 color &＃61; colorRampPalette(c("green", "white", "red"))(100), cellwidth &＃61; 20, cellheight &＃61; 0.2,# 设置单元格大小 fontsize &＃61; 10)

学习思考

——在实战中的问题中(不同分割符号的分割问题)

这是笔者的师兄在做lncRNA数据分析时分享的一个问题&＃xff0c;当时师兄向我分享了先把分割符变成相同的分割符&＃xff0c;然后直接切割这个方法。由于笔者对于separate函数的用法产生兴趣&＃xff0c;help()了一下&＃xff0c;发现原来不同的分割符也可以"一步到位"&＃xff0c;只需要设置sep的参数即可。

所以学习R语言是一个循序渐进的过程&＃xff0c;不可能一蹴而就&＃xff0c;不断学习&＃xff0c;不断提升(共勉&＃xff01;&＃xff01;&＃xff01;)

请多多指教&＃xff01;&＃xff01;&＃xff01;

编辑&＃xff1a;骆栢维

校审&＃xff1a;梁晓杰

相关阅读&＃xff1a;

学习健明老师发布的R语言练习题的学习笔记(一)

利用ROC曲线寻找最佳cutoff值(连续型变量组成的riskscore)

如何使用X-tile软件寻找最佳cutoff值

m6A评分可作为胃癌的预后标志物并预测免疫治疗的效果&＃xff0c;IF&＃61;10.679&＃xff0c;实时IF&＃61;15.83

为何predict()函数计算的Riskscore不等于基因的表达量与其系数的乘积的加权呢?